Curso II – SISSA

Funciones y cómputo en paralelo

Objetivo y contenido

Objetivo y contenido

Objetivo: aprender cómo generar tus propias funciones y realizar cálculos paralelos.

  • Escritura y uso de funciones en R
  • Extracción de precipitación media anual por país
  • Cálculo en paralelo en R
    • Paralelización de procesos

Escritura y uso de funciones en R

Escritura de funciones en R

R permite al usuario crear funciones. El modo de estas expresiones es function, se almacenan en una forma especial y se pueden utilizar en expresiones posteriores.

Una función se define de la siguiente manera:

nombre_funcion <- function(arg1, arg2, ...) expresion

nombre_funcion nombre de la función definida por el usuario.

argn son los diferentes argumentos de la función.

expresion es el proceso que seguirá la función.

Escritura de funciones en R

La función se aplica de la siguiente manera:

nombre_funcion(arg1, arg2, ...)

Escritura de funciones en R

Como ejemplo, generemos una función que cuente los caracteres en una cadena de texto. Primero comenzaremos generando el esqueleto de la función y agregando un parámetro llamado word:

count_char <- function(word){
  
}

Escritura de funciones en R

Luego, aplicamos la función nchar.

count_char <- function(word){
  
  n <- nchar(word)
  
}

Escritura de funciones en R

Finalmente, podemos incluir un mensaje con el resultado. Observa que utilizaremos la función return para indicar explícitamente la salida de nuestra función.

count_char <- function(word){
  
  n <- nchar(word)
  n <- paste("La palabra", word, "tiene", n, "caracteres!")
  return(n)
}

Escritura de funciones en R

Ahora apliquemos la función:

count_char(word = "Oscar Baez")
## [1] "La palabra Oscar Baez tiene 10 caracteres!"
count_char(word = "Análisis espacial")
## [1] "La palabra Análisis espacial tiene 17 caracteres!"
count_char(word = 10:14)
## [1] "La palabra 10 tiene 2 caracteres!" "La palabra 11 tiene 2 caracteres!"
## [3] "La palabra 12 tiene 2 caracteres!" "La palabra 13 tiene 2 caracteres!"
## [5] "La palabra 14 tiene 2 caracteres!"

Writing Functions in R

En el último ejemplo, proporcionamos un vector numérico al parámetro word. Podríamos controlar la entrada del parámetro. Por ejemplo:

count_char <- function(word){
  
  if(!is.character(word))
    stop("¡La entrada no es un caracter!")
  
  n <- nchar(word)
  n <- paste("La palabra", word, "tiene", n, "caracteres!")
  return(n)
}

La función stop interrumpirá la ejecución de la función si no se cumple la condición.

Escritura de funciones en R

count_char("¡Hola chic@s!")
## [1] "La palabra ¡Hola chic@s! tiene 13 caracteres!"
count_char(1528)
## Error in count_char(1528): ¡La entrada no es un caracter!
count_char("¡Este es un curso increíble!")
## [1] "La palabra ¡Este es un curso increíble! tiene 28 caracteres!"

Escritura de funciones en R

¿Cómo crearías una función que te devuelva el valor máximo en un vector dado?

Escritura de funciones en R

¿Cómo crearías una función que te devuelva el valor máximo en un vector dado?

get_maximum <- function(v){
  
  mx <- max(v)
  return(mx)
  
}

Si lo aplicamos:

get_maximum(1:10)
## [1] 10
get_maximum(rnorm(100))
## [1] 2.668295

Aplicación de funciones

Ahora que hemos creado nuestra propia función, ¿cómo podemos aplicarla a un conjunto de datos raster para obtener el valor máximo por celda de grilla durante un período definido? Para este propósito, podemos utilizar la función app del paquete terra.

library(terra)

Aplicación de funciones

Para este propósito, utilizaremos datos mensuales CHIRPSv2 para el año 2000.

chirps_path  <- "C:/User/Data/CHIRPS_monthly"
chirps_files <- list.files(chirps_path, full.names = TRUE)
chirps       <- rast(chirps_files)

Aplicación de funciones

print(chirps)
## class       : SpatRaster 
## dimensions  : 2000, 7200, 12  (nrow, ncol, nlyr)
## resolution  : 0.05, 0.05  (x, y)
## extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : chirps_monthly2000-01.tif  
##               chirps_monthly2000-02.tif  
##               chirps_monthly2000-03.tif  
##               ... and 9 more source(s)
## names       : chirp~00-01, chirp~00-02, chirp~00-03, chirp~00-04, chirp~00-05, chirp~00-06, ... 
## min values  :       0.000,       0.000,       0.000,       0.000,       0.000,       0.000, ... 
## max values  :    1251.456,    1467.558,    1222.492,    1224.301,    1322.237,    1584.747, ...

Aplicación de funciones

Ahora podemos aplicar la función app para obtener el valor máximo en cada celda.

max_values <- app(chirps, get_maximum)
plot(max_values)

Extracción de precipitación media anual por país

Extracción de P media anual

Calcularemos la precipitación media anual utilizando 12 archivos mensuales de CHIRPSv2 sobre los países de la cuenca del Nilo en este ejemplo:

shape_path <- "c:/User/Countries/Nile_Basin_Countries_GAUL2014_2.shp" 
countries  <-vect(shape_path)

Extracción de P media anual

print(countries)
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 15, 9  (geometries, attributes)
##  extent      : 12.19545, 47.98618, -13.45568, 31.66734  (xmin, xmax, ymin, ymax)
##  source      : Nile_Basin_Countries_GAUL2014_2.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :       STATUS DISP_AREA ADM0_CODE   ADM0_NAME STR0_YEAR EXP0_YEAR
##  type        :        <chr>     <chr>     <int>       <chr>     <int>     <int>
##  values      : Member State        NO       205      Rwanda      1000      3000
##                Member State        NO       133       Kenya      1000      3000
##                Member State        NO        74 South Sudan      2011      3000
##  Shape_Leng Shape_Area  Name_label
##       <num>      <num>       <chr>
##       8.127      2.064      Rwanda
##       48.55      47.35       Kenya
##       46.91       51.6 South Sudan

Extracción de P media anual

Por lo tanto, queremos crear una función que:

    1. Tome doce rasters globales de precipitación mensual y genere un stack.
    1. Tome un shapefile con solo un poligono.
    1. Devuelva la precipitación media anual sobre el poligono seleccionado.

Una vez que la función esté lista, se puede aplicar a las diferentes entidades del shapefile.

Extracción de P media anual

Primero, podemos preparar los datos que se utilizarán como entrada en la función. Guardaremos las rutas de los rasters en un objeto llamado raster_paths y el shapefile en un objeto llamado shape.

raster_paths <- "C:/Users/Data/CHIRPS_monthly/"
raster_paths <- list.files(raster_paths, full.names = TRUE)
shape        <- countries[1,]

Extracción de P media anual

print(shape)
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 9  (geometries, attributes)
##  extent      : 28.86187, 30.89965, -2.839881, -1.047913  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :       STATUS DISP_AREA ADM0_CODE ADM0_NAME STR0_YEAR EXP0_YEAR
##  type        :        <chr>     <chr>     <int>     <chr>     <int>     <int>
##  values      : Member State        NO       205    Rwanda      1000      3000
##  Shape_Leng Shape_Area Name_label
##       <num>      <num>      <chr>
##       8.127      2.064     Rwanda

Extracción de P media anual

Luego, podemos crear la función. La llamaremos get_map para la precipitación media anual.

get_map <- function(){
  
}

Extracción de P media anual

Esta función requiere dos entradas: las 12 rutas de los rasters globales de precipitación (uno para cada mes del año) y el shapefile. Podemos proporcionar las rutas de los objetos y leerlos dentro de la función.

get_map <- function(raster_paths, shape){
  
}

Extracción de P media anual

Al final, asignaremos estos objetos a la función, es decir, raster_paths = raster_paths y shape = shape. El primer paso sería importar los datos raster utilizando la función rast del paquete terra.

get_map <- function(raster_paths, shape){
  
  # leyendo los datos
  p_rast <- terra::rast(raster_paths)
  
}

Extracción de P media anual

El segundo paso sería recortar y enmascarar la pila de rasters a la extensión del shapefile.

get_map <- function(raster_paths, shape){
  
  # leyendo los datos
  p_rast <- terra::rast(raster_paths)
  
  # recortando y enmascarando el raster stack
  p_rast <- terra::crop(p_rast, shape, snap = "out")
  p_rast <- terra::mask(p_rast, shape)
  
}

Extracción de P media anual

Finalmente, podemos sumar las 12 capas de precipitación y devolver el resultado.

get_map <- function(raster_paths, shape){
  
  # leyendo los datos
  p_rast <- terra::rast(raster_paths)
  
 # recortando y enmascarando el raster stack
  p_rast <- terra::crop(p_rast, shape, snap = "out")
  p_rast <- terra::mask(p_rast, shape)
  
  # sumando las capas de P
  p_rast <- sum(p_rast)
  
  return(p_rast)
  
}

Extracción de P media anual

Una vez que tenemos la función lista, podemos aplicarla:

rwanda_p <- get_map(raster_paths = raster_paths, shape = shape)
plot(rwanda_p)

Cómputo en paralelo en R

Cómputo en paralelo

El procesamiento de grandes cantidades de datos espaciales puede llevar mucho tiempo. Sin embargo, la mayoría del código de R se ejecuta en un solo procesador. La mayor parte del tiempo esto no es un problema, pero a veces estod procesos pueden:

  • Requerir bastante tiempo

  • Consumir memoria

  • Tomar bastante tiempo al leer o escribir archivos

  • Tomar bastante tiempo al transferir datos

Cómputo en paralelo

Las computadoras tradicionales tienen un solo CPU, que a su vez puede contener múltiples núcleos. Estos procesadores y núcleos pueden realizar cálculos (ten en cuenta que las computadoras modernas pueden tener múltiples procesadores).

Cómputo en paralelo

La siguiente función genera una lista de números primos hasta un número determinado. Esta función fue tomada de stackoverflow.

prime_numbers <- function(n){
   n <- as.integer(n)
   if(n > 1e6) stop("n too large")
   primes     <- rep(TRUE, n)
   primes[1]  <- FALSE
   last.prime <- 2L
   for(i in last.prime:floor(sqrt(n)))
   {
      primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE
      last.prime <- last.prime + min(which(primes[(last.prime+1):n]))
   }
   which(primes)
}

prime_numbers(100)
##  [1]  2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97

Cómputo en paralelo

Utilizaremos la función Sys.time para evaluar el tiempo necesario que lleva esta función para obtener números primos para una secuencia que comienza en 10 y termina en 15,000.

sequence <- 10:15000

# creating an empty vetor
res <- c()

system.time({
  for(i in sequence){
    res[[i]] <- prime_numbers(i)
  }
})
##    user  system elapsed 
##  36.282   0.056  36.351

Cómputo en paralelo

Para paralelizar el código, utilizaremos el paquete doParallel. Este paquete proporciona funciones para la ejecución paralela de código R en máquinas con múltiples núcleos o procesadores o múltiples computadoras.

library(doParallel)

Tenemos que indicar el número de núcleos que se utilizarán. Para ello, podemos usar la función detectCores.

detectCores()
## [1] 12

Cómputo en paralelo

No se recomienda utilizar todos los núcleos para que tu computadora no se bloquee. Por lo tanto:

cores <- detectCores() - 1

Después, tenemos que registrar el backend paralelo. Para ello, utilizaremos la función registerDoParallel.

registerDoParallel(cores = cores)

Cómputo en paralelo

Finalmente, el proceso iterativo del for for será diferente. En lugar de for, utilizaremos foreach; en lugar de in, utilizaremos =. A continuación, podemos indicar explícitamente qué paquetes se utilizarán con el parámetro .packages (repasaremos esto en el ejemplo). Finalmente, el operador %dopar% debe colocarse antes de abrir el foreach.

La función stopImplicitCluster se puede utilizar en documentos y otros lugares donde es importante cerrar explícitamente el clúster creado implícitamente.

Cómputo en paralelo

sequence <- 10:15000

# Parallel 
cores <- detectCores() - 1
registerDoParallel(cores = cores)

# creando un vector vacío
res <- c()

system.time({
  foreach(i = sequence) %dopar% {
      res[[i]] <- prime_numbers(i)
  }
})
##    user  system elapsed 
## 131.530   5.326  14.328
stopImplicitCluster()

Paralelización de procesos

Paralelización de procesos

Retomando el ejemplo de extracción de precipitación media anual por país, obtengamos la precipitación media anual para todos los países de la cuenca del Nilo.

plot(countries, axes = TRUE)

Paralelización de procesos

Al observar los datos, podemos ver que hay áreas en disputa en el shapefile. Como solo queremos los países, podemos excluirlas.

countries <- countries[which(countries$DISP_AREA == "NO"),]
data.frame(countries)
##          STATUS DISP_AREA ADM0_CODE                        ADM0_NAME STR0_YEAR
## 1  Member State        NO       205                           Rwanda      1000
## 2  Member State        NO       133                            Kenya      1000
## 3  Member State        NO        74                      South Sudan      2011
## 4  Member State        NO       257      United Republic of Tanzania      1000
## 5  Member State        NO       253                           Uganda      1000
## 6  Member State        NO        79                         Ethiopia      1000
## 7  Member State        NO        77                          Eritrea      1000
## 8  Member State        NO        43                          Burundi      1000
## 9  Member State        NO        68 Democratic Republic of the Congo      1000
## 10 Member State        NO         6                            Sudan      2011
## 11 Member State        NO     40765                            Egypt      1000
##    EXP0_YEAR Shape_Leng Shape_Area                       Name_label
## 1       3000   8.127003   2.063852                           Rwanda
## 2       3000  48.549430  47.345770                            Kenya
## 3       3000  46.905431  51.599166                      South Sudan
## 4       3000  82.950601  76.860266      United Republic of Tanzania
## 5       3000  23.244416  19.629674                           Uganda
## 6       3000  50.380131  92.869258                         Ethiopia
## 7       3000  50.005783  10.167958                          Eritrea
## 8       3000   9.118459   2.185652                          Burundi
## 9       3000  98.951177 189.998107 Democratic Republic of the Congo
## 10      3000  81.910242 155.888802                            Sudan
## 11      3000  61.251157  89.079113                            Egypt

Paralelización de procesos

Ahora, podemos aplicar la función get_map para cada país. Dado que este proceso seguramente llevará mucho tiempo, podemos paralelizarlo. Primero, aquí está la función get_map.

get_map <- function(raster_paths, shape){
  
  # leyendo los datos
  p_rast <- terra::rast(raster_paths)
  
  # recortando y enmascarando el conjunto de datos ráster
  p_rast <- terra::crop(p_rast, shape, snap = "out")
  p_rast <- terra::mask(p_rast, shape)
  
  # sumando las capas de precipitación
  p_rast <- sum(p_rast)
  
  return(p_rast)
  
}

Paralelización de procesos

Primero, almacenamos las rutas de los archivos ráster en un objeto como hicimos antes. Luego, podemos establecer la ruta de salida donde queremos almacenar los archivos resultantes.

Paralelización de procesos

Construimos el bucle foreach de la siguiente manera:

  • Establecemos el número de núcleos con detectCores
  • Registramos el clúster con registerDoParallel
  • Iniciamos el proceso iterativo foreach
  • Indicamos los paquetes que se utilizarán. En este caso, terra
  • Incluimos el operador paralelo %dopar%
  • Desarrollamos el proceso y exportamos los datos ráster
  • Cerramos el ciclo
  • Cerramos el clúster paralelo si es necesario con stopImplicitCluster

Paralelización de procesos

cores <- detectCores() - 1
registerDoParallel(cores = cores)

foreach(i = 1:nrow(countries), .packages = "terra") %dopar% {
  
  r <- get_map(raster_paths = raster_paths, shape = countries[i,])
  
  n <- countries[i,]$Name_label
  n <- paste0(output, "/", n, ".tif")
  
  writeRaster(r, n, overwrite = TRUE)
  
}
stopImplicitCluster()

Paralelización de procesos

Después del proceso, deberíamos tener los rásteres en la carpeta de salida.

¡Gracias por su atención!